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ABSTRACT 

We present a new semi-analytic model for dynamical friction based on Chandrasekhar’s 
formalism. The key novelty is the introduction of physically motivated, radially vary¬ 
ing, maximum and minimum impact parameters. With these, our model gives an 
excellent match to full N-body simulations for isotropic background density distri¬ 
butions, both cuspy and shallow, without any fine-tuning of the model parameters. 
In particular, we are able to reproduce the dramatic core-stalling effect that occurs 
in shallow/constant density cores, for the first time. This gives us new physical in¬ 
sight into the core-stalling phenomenon. We show that core stalling occurs in the 
limit in which the product of the Coulomb logarithm and the local fraction of stars 
with velocity lower than the infalling body tends to zero. For cuspy backgrounds, this 
occurs when the infalling mass approaches the enclosed background mass. For cored 
backgrounds, it occurs at larger distances from the centre, due to a combination of a 
rapidly increasing minimum impact parameter and a lack of slow moving stars in the 
core. This demonstrates that the physics of core-stalling is likely the same for both 
massive infalling objects and low-mass objects moving in shallow density backgrounds. 
We implement our prescription for dynamical friction in the direct summation code 
NBODY6 as an analytic correction for stars that remain within the Roche volume 
of the infalling object. This approach is computationally efficient, since only stars in 
the inspiralling system need to be evolved with direct summation. Our method can 
be applied to study a variety of astrophysical systems, including young star clusters 
orbiting near the Galactic Centre; globular clusters moving within the Galaxy; and 
dwarf galaxies orbiting within dark matter halos. 

Key words: Galaxy: kinematics and dynamics - Galaxies: star clusters - methods: 
numerical. 


1 INTRODUCTION 


In a seminal work, Chandrasekhar (19431 showed that mas¬ 


sive objects moving through an infinite homogeneous stellar 
medium will experience a drag force that he called ‘dynam¬ 
ical friction’ (see Binney & Tremaine 2008 and . This 


frictional force is likely responsible for galactic mergers (|Gan| 


et al. 2010 Peirani et al. 20101, the accretion of satellites 
and globular clusters onto the Galaxy |Gan et al.| 2010[ 


|Arca-Sedda fc Capuzzo-Dolcetta||2014[ ), and even (in part) 
the coalescence of supermassive black holes (Begelman et al. 


1980). 


Numerical simulations using both direct summation 
codes and collisionless codes have shown that Chan¬ 
drasekhar’s formula (equation]^ works remarkably well, de¬ 
spite the fact that it likely misses important physics like res- 
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onant interactions between the infalling body and the back¬ 
ground (Tremaine fc Weinberg|[i984 Inoue|[2009 Weinberg 


19861. Nonetheless, there are some situations in which it has 


been reported to perform poorly. Most notably, in the case 


of a constant density background (Read et al. 2006 Inoue 
200^|Goerdt et al.|2010T ), or when the mass of the infalling 


body approaches the enclosed background mass ( Gualan¬ 
dris fc Merritt|2008 1. The former is perhaps surprising given 


that the original derivation assumes a homogeneous sea of 
background stars. But this failure could owe to the extreme 
resonance of such harmonic potentials (Read et al.|2006 1. 


While dynamical friction can be accurately modelled 
using N-body simulations, it can often be prohibitively ex¬ 
pensive. To achieve accurate results, the background ’’stars” 
(we shall call them ’’stars”, though these could be any dis¬ 
tribution of gravitating masses, e.g. dark matter particles) 
must be substantially less massive than the infalling body. 
Otherwise, the body will simply be stochastically buffeted 
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around, experiencing little friction ( Baumgardt et al.|2006 |. 
If we wish to self consistently model the internal dynamics 
of a globular cluster, for example, falling onto the Galaxy 
the numerical requirements rapidly become extreme. For 
example a globular cluster of 10®M© (0(10®) particles), 
would need background particles to accurately model 

the inspiral in a massive host. If a lower mass resolution is 
used for the background, and gravitational softening is em¬ 
ployed to keep the relaxation time the same, the dynamical 
friction force is under-predicted (see ^5.2[ ). For this reason, 
semi-analytic models of dynamical friction have become in¬ 
valuable (Hashimoto et al. 2003| Just & Penarrubia 2005 


models computed in NBODY6 and the tree-code GADGET 


Just et al.|2011 Arca-Sedda et al. 2015 I. These significantly 

speed up the simulations since only the internal dynam¬ 
ics of the satellite needs to be integrated self-consistently, 
and the effects of a particular background distribution are 
modelled in an analytic way. Such models have been well- 
tested for point mass satellites in steep power-law density 
backgrounds, giving a good match to full N-body simula¬ 
tions (Just et al. 20111. However, in shallow or constant 


density backgrounds, dynamical friction stalls (Read et al. 


2006) - an effect that has so far not been captured by 


semi-analytic models. This has led some authors to move 
away from semi-analytic modelling, towards particle-mesh 


codes calibrated by direct Y-body and tree-codes (Spin- 
nato et 5I||2003), mixed collisional/collisionless methodolo¬ 


gies ( Fujii et al.||2006| 2007[ [2009 1, and accurate advanced 
tree-codes | Dehnen|2014[ ). While these are exciting new ap¬ 
proaches, they remain computationally expensive. This begs 
the question of whether the semi-analytic models cannot be 
improved. Such improvements would not only open up new 
classes of astrophysically interesting problems, but also shed 
more light on the interesting physics of dynamical friction 
and core stalling. This is the goal of this present work. 

In this paper, we introduce a new semi-analytic model 
for dynamical friction based on the familiar Chandrasekhar 
formalism (equation [^. Our key novelty is that we present 
a new physically motivated, and radially varying minimum 
impact parameter (6min), which when combined with the 
maximum impact parameter (6, 


from Just & Penarrubia 


120051, gives our semi-analytic model a remarkable match to 


a large range of full N-body simulations, without any need 
for fine tuning of the model parameters. We introduce an 
ansatz that dynamical friction ceases when bmin > 5max, as 
there are no valid impact parameters for encounters that 
would contribute to the frictional force (see section section 
6.11. We describe and test a method of applying equation 
to a cluster of particles, and implement this model in the di¬ 
rect summation code NBODY6 ( Nitadori fc Aarseth|201^ . 
This opens up a wide array of astrophysically interesting 
problems, for example following the internal dynamics of 
collisional star clusters while simultaneously capturing their 
orbital decay due to dynamical friction in the inner galaxy. 

The paper is organised as follows. In section il we de¬ 
scribe the physical motivation for our choice of impact pa¬ 
rameters and the treatment of fcmax in cored profiles. Sec¬ 
tion describes how Ghandrasekhar’s formula is applied 
to satellites comprised of a cluster of particles and its im¬ 
plementation in NBODY6. Section describes the galaxy 
and cluster models used throughout the paper. In section il 
we compare the semi-analytic results of our modified ver¬ 
sion of NBODY6 (hereafter NBODY6df) with full Y-body 


(Springel et al. [2001 1 . Finally we compare with previous 


work in section ^ and conclude our results in section 0 


2 THEORY 


Chandrasekhar’s dynamical friction formula for a satellite of 
mass Ms is often written as (Binney fc Tremaine||20d8 ): 


^ = -47rG^Msplog(A)/(u.)^, 
at Vs 


( 1 ) 


where log(A) is the Coulomb logarithm equal to 
log (6inax/5min), vg is the Satellite velocity, us = [usl and p 
is the local background density. If we assume a Maxwellian 
distribution of velocities: 

9 Y 

/(u*) = erf(X)-^exp(-A:^), (2) 

where X = vs/-\/2a. 

Equation!^ is derived under the assumption of an infi¬ 
nite homogeneous background. Despite this assumption, the 
Coulomb logarithm, log(A), takes into account the finite size 
of a real system through the maximum and minimum im¬ 
pact parameters. Typically 6max is of the order of the size of 
the host system, and bmin is defined as the impact parameter 
for a 90 degree deflection. 

From a theoretical standpoint, it is surprising that 
Chandrasekhar’s formula is so successful at reproducing the 
effects of dynamical friction. In real systems, dynamical fric¬ 
tion almost certainly results from discrete resonances with 


background stars (Tremaine & Weinberg 1984 Weinberg 


1986). Chandrasekhar’s formula does not model these res¬ 


onances, however it likely works because in most systems 
when an infalling body migrates from one radius to the next 
it experiences a whole new set of resonances. [Tremaine fc| 
Weinberg (1984 1 show that if the resonances are assumed 


to form a continuum, Chandrasekhar’s approximation is re¬ 
produced. In this way the infalling body behaves similarly 
to a massive body moving through an infinite homogeneous 
medium that encounters each background star only once. If 
we construct a special system, however, where by moving 
from one radius to the next we do not encounter new res¬ 
onances then Chandrasekhar’s formula has been known to 
fail. An example of this is given by the harmonic core (where 
the background density is constant: p(r) ~ po), which ex¬ 
hibits a short burst of super-Chandrasekhar dynamical fric¬ 


tion followed by rapid stalling of the cluster orbit (Read 
et al.|2006 Inoue|2(309 Goerdt et al.|20l0 |. 


2.1 Maximum impact parameter 

Although a constant 6max has been traditionally used to esti¬ 


mate inspiral timescales, Hashimoto et al. (20031 computed 


semi-analytic approximations of Y-body models and found 
that a spatially dependent bmax better reproduces simula¬ 
tion results. The physical motivation for this comes from 
the local approximation under which equation is derived, 
which assumes the density distribution is constant up to 
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&max- Therefore 6 max should be a local property of the sys¬ 
tem. The authors took &max to be the distance from the 


of magnitude estimate, as the density of particles with im¬ 
pact parameters larger than i?g is low compared to the local 
density distribution around the subject cluster. However for 
sufficiently cuspy profiles (7 > 1 , where 7 is the asymptotic 
slope of the distribution) this approach will typically over 
estimate the dynamical friction effect near the centre of the 
background distribution. The slope of the density distribu¬ 
tion is difficult to account for, however Just et al. (20111 
show that: 


f^max — 


P(^g) 

dp I 
'dr \ 


(3) 


is a better maximum cutoff to compensate for the cuspy den¬ 
sity profile (i.e. the local density over the local density gra¬ 
dient). This impact parameter gives a length scale for which 
the density is approximately constant, giving a more accu¬ 
rate representation of the local approximation (see also |Just| 
fc Penarrubia|[2b05 1. This makes intuitive sense if one con¬ 
siders the two extreme cases of density distribution. If p is 
a constant over all space, fomax —>■ cio, and the force logarith¬ 
mically diverges, as in Chandrasekhar’s original derivation. 
On the other hand if the distribution is infinitely cuspy (i.e. 
p{r) ~ J'“°“), femax —> 0 , and the satellite effectively orbits a 
point mass in a Keplerian orbit with no decay. 

From equation [j} |Just et al. (20111 approximate that 
&max ~ Rg/^j however we keep the full expression in this 
paper so that the denominator can vary locally. We as¬ 
sume that the background density distribution is given by 
a Dehnen model ( |Dehnen[]1993p . For this model the density 
and its derivative are both analytic, and Equation can be 
expressed as: 


f^max — 


07 -I- iRg 


(4) 


which indeed reduces to 6 max = RgM for Rg « a. Note 
that an attractive feature of equation is that it is well- 
behaved in the limit 7 —>■ 0 , tending to a constant femax = 
a/4. This is not the case if we use instead femax = Rg/j 
for which 6max —>■ 00 . This led Just et al. (20111 to adopt 


= Rg for 7 < 1. It would seem, then, that our equation 
l^is superior in this regard. However, the finite 7 —> 0 limit 
is peculiar to the assumed background Dehnen profile. It is 
straightforward to show that split power law profiles that 
transition from the inner to the outer slope more steeply 
than the Dehnen profile will a produce divergent fomax in the 
limit 7 —> 0 , if we assume the ansatz 6max = p{Rg)/^PiRg)- 


For this reason, in this paper we follow Just et al. (20111 
and assume 6max = Rg for 7 < 1 . 


2.2 Minimum Impact Parameter 

The minimum impact parameter (i.e. the impact parameter 
corresponding to a 90 degree deflection) of extended objects 
is roughly of the order of the half mass radius of the object 
( Binney fc Tremaine||200'8 1 . Hashimoto et al. (20031 found, 
for infalling satellites of a Plummer density profile, that bmin. 
is well approximated by 1.4es, where ts is the Plummer scale 



Figure 1. Schematic of the velocity distribution function. If 
iS O', then [Just et al. ] pMTj l’s Wtyp significantly overesti¬ 
mates the maximum impact velocity, as, for isotropic distribu¬ 
tions, background stars moving faster than the satellite (the 
shaded area) are not expected to induce dynamical friction. The 
schematic assumes a Maxwellian velocity distribution and is nor¬ 
malised to units of (7. 


radius, a. In terms of the half-mass radius this corresponds 
to foniin = (1.4/1.3)rhm. 

It should be noted that even though jHashimoto et al.| 
( [2003] ) ’s 6 min was fit for a Plummer sphere, it is a reason¬ 
able approximation for other stellar distributions. [Just fc| 
Penarrubia (20051 show in their discussion about 6 min that 


the minimum impact parameter for Plummer, King and sin¬ 
gular isothermal sphere models are very similar in terms of 
the half mass radius. Similarly to Just & Penarrubia (20051, 
we take fomin of extended to objects to be rhm instead of 
Hashimoto et al. (2003l’s ~ 1.07rhm, to keep our formal¬ 


ism physically motivated rather than calibrated by V-body 
models. 

At any epoch, we take fcmin to be the maximum of rhm 
and the minimum impact parameter of a point mass, which 
is typically taken to be GM/vtyp^ (where M and Vtyp are the 
bound mass and ’’typical” velocity for an encounter, respec¬ 
tively) ( [Binney Sz Tremaine|2008[). In general, uty p (and thus 
5min) are poorly constrained. Just & Penarrubia ( |2005| ) take 
Wtyp^ = 2 (t^ -|- us^. If one considers that b will be minimised 
for the highest velocity encounter, this seems a reasonable 
choice at first glance. However, with this formulation Utyp is 
the maximum relative velocity of any encounter (i.e. the tail 
of the Maxwellian velocity distribution), not the maximum 
velocity of encounters that actually contribute to dynamical 
friction. With this prescription when ~ (which can 
be the case at the apocentre of eccentric orbits, or through¬ 
out the entire core in shallow profiles) the minimum im¬ 
pact parameter will be largely underestimated. This is be¬ 
cause the relative velocity of the encounter is larger than the 
velocity of the satellite, and (for isotropic background dis¬ 
tributions) only particles moving slower than the satellite 


are assumed to contribute to dynamical friction (Binney & 
Tremaine|2008 1 . The Maxwellian velocity distribution func¬ 


tion already selects only the stars that have velocities slower 
than the satellite to contribute to the friction, and as such 
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only the relative velocities of these stars need to be taken 
into account, for consistency. For these reasons, we propose 
that Vtyp should really be interpreted as vs, the velocity of 
the infalling body. This is the maximum velocity encounter 
that can contribute to dynamical friction and will therefore 
give rise to the smallest impact parameter. As such we re¬ 
define Vtyp^ = vs^- Fig. displays this schematically. The 
implications of using this corrected minimum impact pa¬ 
rameter over Just et ^ ( |2011 1 ’s prescription are discussed 
in section l5.2.2l 


2.3 The Coulomb logarithm and core stalling 

These prescriptions for bmax and femin give us the following 
functional form for the Coulomb logarithm: 


(5) 

Equation shows that our prescription for log A is a func¬ 
tion of the radial distance to the centre of the background 
potential, the slope of the background distribution and the 
half mass radius of the cluster. If during inspiral femin ^ &max, 
the dynamical friction term is set to zero. This ansatz is rea¬ 
sonable since this means there are no particles available to 
scatter off the satellite in a way that would reduce its orbital 
energy We now show that this ansatz is equivalent to the 
well known result that friction ceases if the satellite mass 


e|2008D , 


approaches the enclosed mass of the background (Binney & 
|Trer 


Vtyp 


2 gmarA gmarA 


Rg 

GM 

Vtyp^ 

M 

MARA’ 


M 


Mg{Rg 


(6a) 

(6b) 

(6c) 


where Mg{Rg) is the galaxy mass enclosed at Rg. 

Stalling occurs at this scale because perturbations from 
individual stars dominate over the mean field effects, making 


dynamical friction less efficient (Gualandris & Merritt 20081. 


For the case of a Dehnen model an approximate analytic 
equation for the stalling radius can be derived. Equating the 
argument of the Coulomb logarithm to unity and assuming 
a circular orbit: 


^max Rg{Rg + a)/4:Rg 

b~ ^ GMs/vA ’ 

^2 GMg{Rg) 

V s = -^-, 

Kg 

^max _ Rg{Rg + a)/4:Rg 

^min MsRg/Mg{RA ’ 

^max _ ^ Rg + / Mg{Rg) \ _ 

femin ^\ARg)\ Ms J ^ 


Recalling the formula for Mg{Rg) ( Dehnen]|1993 1: 


Mg{RA = Mg 


Rg 


Rg CL 


3-7 


(7a) 

(7b) 

(7c) 

(7d) 

(7e) 


and inserting this into equation |7d| and rearranging: 


Ms 

f Rg y-^Rg + a 

(7f) 

Mg 

\Rg + a J 4Rg 

Ms 

A *' 

(7g) 

Mg 

(r + aA (07 + 4r) 

If we take the limit of r < a: 


Ms 

A-f 

(7h) 

Mg 

A — 1 J- 07 

Therefore: 



Rs = 

( Ms , 2-7 , 

(7i) 


where Rs is the stalling radius of the satellite. Note that 
this is the theoretical stalling radius for a point particle. If 
the particle loses mass, Ms/Mq will shrink and the cluster 
can potentially reach further in, but of course the timescale 
for inspiral will be longer. We will show in section |5.2.2| 
that this shrinking log (A) captures the unique core stalling 
in shallow profiles when coupled with /(v*). For a profile 
with an intrinsically fiat core, stalling occurs even farther 
out, this phenomenon is explained in section [6T| 


2.4 Velocity Dispersion 


The fraction of background stars moving slower than the 
satellite (equation]^ is obtained from the underlying den¬ 
sity distribution. Given a particular analytic density distri¬ 
bution, the velocity dispersion as a function of Rg can be 
derived from the Jeans equation. For a Dehnen model cr(Rg) 
is analytic if dy is an integer, and in the current implementa¬ 
tion of NBODY6df a selection of analytic results have been 
implemented for various values of gamma. Once derived this 
allows for a quick analytical calculation in the code (see Ap¬ 
pendix for the full derivation). To use a non-integer value of 
47 one would need to implement a numerical solver of the 
Jeans equation in the code. However, for the sake of speed, 
we suggest instead to look for a degenerate model by modify¬ 
ing the scale radius, a, and mass, Mq, of the Dehnen model 
so that an integer 47 may be used. We chose to use Dehnen 
models in the initial implementation due to their versatility 
for modelling spherical systems. If one would like to imple¬ 
ment a different density distribution, we urge authors to 
take great care with the definition of 6max = po/Ap{r)o, 
however this impact parameter is calculable for any density 


distribution with a continuous derivative (see also Just & 
Penarrubia (20051; Just et al. ( 2011| ). 


3 IMPLEMENTATION 


We implemented equation as an external analytical accel¬ 
eration in NBODY6df. Dynamical friction is applied on the 
regular integration step, which is computed in parallel on 
the GPU. 


Fujii et al. (20061 fit semi-analytical models to V-body 


simulations of dwarf galaxies experiencing dynamical fric¬ 
tion in larger parent galaxies. It is shown that the simula¬ 
tions undergo enhanced dynamical friction as compared with 
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Chandrasekhar’s formula due to two effects. The first is di¬ 
rect gravitational interactions with escaped particles. This is 
included in NBODYGdf by integrating tidally stripped ma¬ 
terial self consistently. The second is the indirect effect of 
material that is energetically unbound but remains close to 
the cluster, enhancing the gravitational wake in the back¬ 
ground medium. In an attempt to replicate this effect, the 
mass Ms in equation is taken to be the total mass of 
the particles contained inside the Roche volume, i.e. the 
particles for which Fciuster > Ttidai, where Fduster is the 
magnitude of the gravitational force on the star due to the 
A^-body particles and Ttidai is the tidal force experienced by 
the star. Ftidai is defined as IT’g —i^o|, where is the force 
experienced by the star due to the background distribution, 
and Fo is the similar force felt by the density centre of the 
cluster. This procedure ensures correct calculation of Roche- 
volume membership and is scale independent, requiring no 
approximate tidal radius. The force exerted by the back¬ 
ground distribution is calculated using existing NBODY 6 
routines. If the user wishes to change the background poten¬ 
tial they can simply replace the potential in NBODYGdf in 
the same manner they would in NBODY 6 , however we again 
stress that the velocity dispersion calculation and maximum 
impact parameter need to be updated too, as discussed in 
the previous section. We will refer to particles inside the 
Roche volume as “bound” for ease of use, even though this 
includes potential escapers. Particles experience no dynam¬ 
ical friction whilst unbound, but will feel dynamical friction 
once again if they re-enter the Roche volume. NBODY 6 in¬ 
cludes the regularisation of binaries and close encounters 
( Mikkola fc Aarseth||1998 1 in which the system is replaced 
by a centre of mass (CoM) particle during the regular step. 
The regularised system is considered under the dynamical 
friction regime if its CoM particle is bound. The dynamical 
friction force is then applied directly to the CoM particle, 
and the differential force on each particle is handled by the 
KS regularisation scheme. 

Collecting all of the scalar terms in equation [fallows it 
to be rewritten as: 


®df — C’fric('Us/u s). 


( 8 ) 


where: 


C^fric — 47rMRoche log(A)p(Rg) 


erf(X) 


2X 

7 ^ 


exp(—Y^) 


(9) 

The cluster velocity, vq, is taken to be the average velocity 
of the particles in the cluster core with respect to the galac- 
tocentric rest frame. Equation is calculated every time 
NBODY 6 df adjusts important parameters. This coefficient 
is used in all dynamical friction calculations until the next 
adjust routine. The coefficient varies slowly and the com¬ 
puted value is approximately constant between parameter 
adjustments so long as the adjust time is significantly lower 
than the orbital period. 

When the dynamical friction correction is applied to a 
particle, the change in energy is calculated and added to 
a variable representing the total energy removed from the 
system due to dynamical friction. When NBODYGdf calcu¬ 
lates the total energy of the system we include this term so 
that total energy is conserved, and the energy error in the 
Y-body calculation can be evaluated in the usual way. 


4 SIMULATIONS 
4.1 Initial Conditions 


4-1.1 Cluster 

The clusters in this study are initially Plummer models of 
mass Ms = 10® Mq and half mass radius rhm = 0.1 pc, 
similar to the clusters modelled in Kim & Morris (20031 
and Fujii et al. (20091. The mass of a cluster particle is 
ms = Ms/Ns, where Ns is the number of cluster particles. 


4.1.2 Background 


We adopt single component Dehnen models | Dehnenjl993 1, 
representing the central region of the Galaxy. We use a slope 
7 = 1.5, scale radius a = 8.625 pc and mass Mg = 5.9 x 
10 ^ Mq to represent the density distribution in the central 
few tens of parsecs in the Milky Way. This closely represents 
the observed broken power-law profile obtained by |Genzel| 
et al. (20031 for the central 10 parsecs of the Galaxy. For 


runs with different 7 we use the same parameters as stated 
above. What these profiles represent is arbitrary, and we 
simply keep the same mass and scale radius for ease of use. 

It should be noted that any time-independent analyt¬ 
ical spherical background potential can be included in the 
code, such as the addition of a central SMBH, and a dark 
halo component (although for some models the density and 
velocity dispersion functions may need to be calculated nu¬ 
merically). We have adopted a simple model here to ease 
comparison of the code with full Y-body models with low- 
Y, in which a SMBH may wander signihcantly. Ghoosing a 
single spherical component for the test simulations also gives 
more applicability to larger scale simulations. Models where 
the scale radius and mass of the host are larger (i.e. glob¬ 
ular clusters in a dwarf spheroidal) will behave in a similar 
way, as integration in NBODY 6 is performed in scale inde¬ 
pendent Henon units internally, G = 1 = M = = —iE 

(where M is the total mass, Rv is the virial radius and E is 


the total energy) (Nitadori & Aarseth 20121. A more real 


istic treatment of clusters near the Galactic Centre will be 
investigated in a scientihe context in an upcoming paper. 


4.2 Models 

We compare results of NBODY 6 df with results from fully 
self-consistent NBODY 6 and GADGET runs, where the 
background distribution is granular. NBODY 6 is a direct- 
summation collisional code and as such we use equal particle 
masses for the cluster and background to reduce unrealistic 
scattering. When using the tree-code GADGET however, 
we use a smaller mass for the cluster particles, as softening 
reduces the low-Y scattering effects. We use softening pa¬ 
rameters of 0.025 pc for the cluster particles and 0.1 pc for 
the bulge particles in all GADGET simulations. 

With the maximum initial cluster distance being 10 pc 
and the background i?hm being ~ 15 pc for 7 = 1.5, the 
truncation of the background at large radii will have a neg¬ 
ligible effect on the dynamical friction experienced by the 
cluster. Thus for both the GADGET and NBODY 6 runs 
we truncate the Dehnen potential at 50 pc. The models are 
summarised in Table [T] 
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Run Name 

Code 

Ns 

m*S 

^bg 

m*bg 

Ra. 

7 

Eccentricity 

Section 




(M©) 


(M©) 

(pc) 




dflk 

NBODYOdf 

Ik 

100.0 

- 

- 

10.0 

1.5 

0.0 

5.1.1 


dflke 

NBODYOdf 

Ik 

100.0 

- 

- 

10.0 

1.5 

0.3 


dflkeO.75 

NBODYOdf 

Ik 

100.0 

- 

- 

10.0 

1.5 

0.75 


dflkgl 

NBODYOdf 

Ik 

100.0 

- 

- 

10.0 

1.0 

0.0 


dflkgl.75 

NBODYOdf 

Ik 

100.0 

- 

- 

10.0 

1.75 

0.0 


nblk 

NBODY6 

Ik 

100.0 

464k 

100.0 

10.0 

1.5 

0.0 


nblke 

NBODY6 

Ik 

100.0 

464k 

100.0 

10.0 

1.5 

0.3 


nblkeO.75 

NBODY6 

Ik 

100.0 

464k 

100.0 

10.0 

1.5 

0.75 


nblkgl 

NBODY6 

Ik 

100.0 

429k 

100.0 

10.0 

1.0 

0.0 


nblkgl.75 

NBODY6 

Ik 

100.0 

483k 

100.0 

10.0 

1.75 

0.0 


dfL5 

NBODYOdf 

Ik 

100.0 

- 

- 

5.0 

1.5 

0.0 

5.1.2 


dtL5e 

NBODYOdf 

Ik 

100.0 

- 

- 

5.0 

1.5 

0.3 


dfL5.0e0.75 

NBODTOdf 

Ik 

100.0 

- 

- 

5.0 

1.5 

0.75 


dfL2.5 

NBODYOdf 

Ik 

100.0 

- 

- 

2.5 

1.5 

0.0 


dfL2.5e 

NBODYOdf 

Ik 

100.0 

- 

- 

2.5 

1.5 

0.3 


dfL2.5eO.75 

NBODYOdf 

Ik 

100.0 

- 

- 

2.5 

1.5 

0.75 


nbL5 

NBODY6 

Ik 

100.0 

464k 

100.0 

5.0 

1.5 

0.0 


nbL5e 

NBODY6 

Ik 

100.0 

464k 

100.0 

5.0 

1.5 

0.3 


nbL5.0e0.75 

NBODY6 

Ik 

100.0 

464k 

1000.0 

5.0 

1.5 

0.75 


nbL2.5 

NBODY6 

Ik 

100.0 

464k 

100.0 

2.5 

1.5 

0.0 


nbL2.5e 

NBODY6 

Ik 

100.0 

464k 

100.0 

2.5 

1.5 

0.3 


nbL2.5eO.75 

NBODY6 

Ik 

100.0 

464k 

1000.0 

2.5 

1.5 

0.75 


dfL2k 

NBODYOdf 

2 k 

50.0 

- 

- 

5.0 

1.5 

0.0 


dfL2ke 

NBODYOdf 

2 k 

50.0 

- 

- 

5.0 

1.5 

0.75 


dfL4k 

NBODYOdf 

4k 

25.0 

- 

- 

5.0 

1.5 

0.0 


dfL4ke 

NBODYOdf 

4k 

25.0 

- 

- 

5.0 

1.5 

0.75 


dfL8k 

NBODYOdf 

8 k 

12.5 

- 

- 

5.0 

1.5 

0.0 


dfL8ke 

NBODYOdf 

8 k 

12.5 

- 

- 

5.0 

1.5 

0.75 


nbL2k 

NBODY6 

2 k 

50.0 

464k 

100.0 

5.0 

1.5 

0.0 


nbL2ke 

NBODY6 

2 k 

50.0 

464k 

100.0 

5.0 

1.5 

0.0 


nbL4k 

NBODY6 

4k 

25.0 

464k 

100.0 

5.0 

1.5 

0.0 


nbL4ke 

NBODY6 

4k 

25.0 

464k 

100.0 

5.0 

1.5 

0.75 


nbL8k 

NBODY6 

8 k 

12.5 

464k 

100.0 

5.0 

1.5 

0.0 


nbLSke 

NBODY6 

8 k 

12.5 

464k 

100.0 

5.0 

1.5 

0.75 


dfalO 

NBODYOdf 

lOk 

10.0 

- 

- 

10.0 

1.5 

0.0 

5.2 


dfalOe 

NBODYOdf 

lOk 

10.0 

- 

- 

10.0 

1.5 

0.3 



dfaS 

NBODYOdf 

lOk 

10.0 

- 

- 

5.0 

1.5 

0.0 


dfa2.5 

NBODYOdf 

lOk 

10.0 

- 

- 

2.5 

1.5 

0.0 


gtalO 

GADGET 

lOk 

10.0 

1549k 

30.0 

10.0 

1.5 

0.0 


gtalOe 

GADGET 

lOk 

10.0 

1549k 

30.0 

10.0 

1.5 

0.3 


gtaS 

GADGET 

lOk 

10.0 

1549k 

30.0 

5.0 

1.5 

0.0 


gta2.5 

GADGET 

lOk 

10.0 

1549k 

30.0 

2.5 

1.5 

0.0 


gtagO.O 

GADGET 

lOk 

10.0 

1549k 

30.0 

5.0 

0.0 

0.0 

5.2.2 


gtagO.5 

GADGET 

lOk 

10.0 

1549k 

30.0 

5.0 

0.5 

0.0 


dfgO.O 

NBODYOdf 

lOk 

10.0 

- 

- 

5.0 

0.0 

0.0 


dfgO.5 

NBODYOdf 

lOk 

10.0 

- 

- 

5.0 

0.5 

0.0 


df2k 

NBODYOdf 

2 k 

50.0 

- 

- 

10.0 

1.5 

0.0 

5.3 


df4k 

NBODYOdf 

4k 

25.0 

- 

- 

10.0 

1.5 

0.0 



df8k 

NBODYOdf 

8 k 

12.5 

- 

- 

10.0 

1.5 

0.0 


dfl6k 

NBODYOdf 

16k 

6.25 

- 

- 

10.0 

1.5 

0.0 


df32k 

NBODYOdf 

32k 

3.125 

- 

- 

10.0 

1.5 

0.0 


df64k 

NBODYOdf 

64k 

1.5625 

- 

- 

10.0 

1.5 

0.0 



Table 1. Initial conditions of simulations. Column 1 lists the names of the simulations, where the prefixes: df, nb and gt indicate 
the code used, NBODYhdf, NBODY 6 and GADGET, respectively. This is also stated in column 2. Golumns 3 to 6 display the 
particle numbers and masses for both the cluster and the background, subscripts c and bg respectively. Column 7 lists the initial 
distance of the cluster from the Galactic Centre, all runs start at apocentre. Column 8 states the asymptotic slope used in the 
background Dehnen model. Column 9 shows the initial eccentricity of the cluster. Column 10 displays which chapter each group of 
simulations first appears in. 
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5 RESULTS 

5.1 Comparison with NBODY6 

5.1.1 Orbit Comparison 


Simulations dflk and nblk are compared in Fig. which 
shows the radial position of the cluster with respect to the 
Galactic Centre as a function of time. Fig.|^shows the bound 
mass of the clusters in the different simulations. The agree¬ 
ment between the two models is excellent. After ~ 2 Myrs 
nblk experiences stochastic changes in its orbit due to the 
\ow-N background, this is because the low-N cluster has 
nearly dissolved by this time. Prior to this epoch, when the 
clusters are not close to dissolution, the radial distance trav¬ 
elled by the cluster in the two codes differs by less than a 
few per cent. 

The stochastic changes in the cluster orbit of simula¬ 
tion nblk come from Y-body sampling from the distribu¬ 
tion function. This introduces chaotic effects on both large 
and small scales, compared to the equivalent analytic dis¬ 
tribution ( Arca-Sedda fc Capuzzo-Dolcetta||2014 l. At small 
scales the granularity of the background induces stochastic 
changes in the orbit if N of the background is low (in nblk 
each background particle represents the mass of an ensem¬ 
ble of stars). On large scales the system may deviate from 
spherical symmetry, inducing moderate eccentricity and pre¬ 
cession. These effects accumulate over time and cause the 
eccentricity in nblk once the cluster has almost dissolved. 

Fig-El shows snapshots of the simulations at different 
times. Only after ~ 3 Myrs are the models distinguishable, 
and it can be seen that the structure and distribution of the 
tidally stripped material is well reproduced in NBODYGdf. 

Figs. 0 and show simulations which have the same 
initial conditions as dflk and nblk, but with initial cluster 
eccentricities of 0.3 and 0.75, respectively. The agreement 
is excellent for both eccentricities. For e = 0.75 the simula¬ 
tions diverge near the end as the clusters have lost the ma¬ 
jority of their mass. The agreement is so good because of our 
prescription for femin, which is dependent on both position 
and velocity. At pericentre the cluster moves fastest, giving 



Figure 2. Distance of the cluster with respect to the Galactic 
Centre as a function of time for dflk (blue line) and nblk (red 
dashed line). 



Figure 3. Mass enclosed in the Roche volume as a function of 
time for dflk (blue line) and nblk (red dashed line) 
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Figure 4. Snapshots of the models dflk (blue) and nblk (red) 
at different epochs. For nblk only particles originating from the 
cluster are plotted. 


a smaller bmin and a stronger dynamical friction force, at 
apocentre the opposite is true, decreasing the force. Mean¬ 
while bmax varies across the length of the orbit due to its 
radial dependence. The result is an accurate calculation of 
the force along the entire orbit. The excellent agreement of 
both these models shows that the semi-analytic dynamical 
friction scheme in NBODYGdf can accurately reproduce the 
force experienced for a range of eccentricities. 

Figures and compare models which have the same 
initial conditions as dflk and nblk, but with asymptotic 
slopes of 7 = 1.75 and 7=1 respectively. Both show good 
agreement. The NBODY 6 runs gain some eccentricity from 
the granularity of the low-Y background distributions. This 
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Figure 5. Distance of the cluster with respect to the Galactic 
Centre as a function of time for dflke (blue line) and nblke (red 
dashed line). 



Figure 7. Distance of the cluster with respect to the Galactic 
Centre as a function of time for dflkgl75 (blue line) and nblkgl75 
(red dashed line). 



8.0 0.5 1.0 1.5 2.0 2.5 

t (Myr) 

Figure 6. Distance of the cluster with respect to the Galac¬ 
tic Centre as a function of time for dflkeO.75 (blue line) and 
nblkeO.75 (red dashed line). 

common problem with the low A^-models is addressed in 
sections [ 5 T 21 and 15.21 

5.1.2 Angular Momentum Comparison 

During inspiral the dynamical friction force is coupled 
with the relaxation of the cluster, meaning &min(f?g,us) = 
buiin{Rg,vs,t) fl'iid Mci = Mci(t). Therefore different reali¬ 
sations of low-A'^ simulations can significantly deviate from 
each other by using a different random seed, as the mass loss 
from dynamical ejections is very much a stochastic process 
for low-A^ simulations, where the relaxation time is short. 

Attempting to isolate each effect can give some indica¬ 
tion of how accurate an approximation NBODYGdf is. In the 
limit of negligible dynamical friction NBODY6df is identical 
to NBODY6, and the relaxation timescales will be similar. 

We ran a series of short simulations to try to isolate 
the dynamical friction effect from the relaxation process as 



Figure 8. Distance of the cluster with respect to the Galactic 
Centre as a function of time for dflkgl (blue line) and nblkgl 
(red dashed line). 

much as possible. In Fig. we plot the total angular mo¬ 
mentum of the bound material perpendicular to the orbital 
plane (i.e. LA as a function of time for half an orbit, for 
different initial cluster orbits (see table for initial condi¬ 
tions). Fig. |10| shows the same for eccentricities of 0.75. Over 
a time of only half an orbit the clusters lose no more than 
10 per cent of their mass, and as such the orbital evolution 
is only weakly dependant on relaxation. In Fig. the full 
Ai-body models lose a bit more mass than NBODYOdf at 
pericentre. After further investigation this seems to be due 
to stars near the tidal radius of the cluster being stripped 
more aggressively due to two body scattering with the low-A’ 
background, rather than a deviation from Chandrasekhar’s 
formula. 

We tested a specific case (a=5.0,e=0.0,0.75) in models 
where the cluster is comprised of 2k, 4k and 8k particles, 
whilst keeping the total mass constant. We did not redo the 
entire grid of initial radii and eccentricities as the full N- 
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0.2 0.4 0.6 0.8 

t (Myr) 


Figure 9. Z-component of angular momentum (perpendicular to 
the orbital plane) as a function of time for half a cluster orbit at 
different initial distances 10, 5 and 2.5 pc, and eccentricities of 
0.0 and 0.3. In all cases 7 = 1.5. 


Figure 11. Z-component of angular momentum (perpendicular 
to the orbital plane) as a function of time until cluster loses 10 % 
of its mass. At an initial distance of 5 pc, witch eccentricities of 
0.0 and 0.75. From top to bottom the cluster consists of 2k, 4k 
and 8 k particles. In all cases 7 = 1.5. 


6.0 


dflkeO.75 &L nblkeO.75 



A 

'i/i 

E 

O 3.0 


2.0 


1.0 


dfL5e0.75 & nbL5e0.75 


dfL2.5eO.75 & nbL2.5eO.75 


- NBODY6df 

-- NBODY6 


circular and eccentric orbits. The orbital evolution over 
many orbits can be considered accurate because the dynam¬ 
ical friction coefficient is linear with mass. At any epoch, 
ti, Ms = Ms(ti) and Rg = Rg{ti). If we assume that our 
limited number of models in Fig. indicate that the dy¬ 
namical friction coefficient is initially correct compared to 
full A-body models at any Ms,o, Rg,o, eo, then an entire in¬ 
spiral can be thought of as traversing a grid of these models, 
and as such the dynamical friction coefficient, when decou¬ 
pled from relaxation, can be considered correct. The models 
in Fig. have an asymptotic slope of 7 = 1.5, but similar 
agreement is found for 7 = 1 and 7 = 1.75. 


°tPoo 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 

t (Myr) 

Figure 10. Z-component of angular momentum (perpendicular 
to the orbital plane) as a function of time for half a cluster orbit 
at different initial distances 10, 5 and 2.5 pc, with eccentricity of 
0.75. In all cases 7 = 1.5. 


body models are very numerically expensive. The models 
were run until the cluster in the NBODY6df simulation had 
lost 10% of its mass. Fig. shows these sets of models, 
where very good agreement is found. The higher N mod¬ 
els retain their mass for longer due to slower two-body re¬ 
laxation. In all the NBODY6 simulations, the high ratio of 
cluster particle mass to background particle mass causes a 
few of the less tightly bound cluster particles to be stripped 
very early on in the simulation. Therefore the angular mo¬ 
mentum of the bound particles is systematically lower af¬ 
ter ~ 0.2 Myr in all cases. This is not significant, as these 
stripped particles contribute very little to the total mass of 
the cluster, and so the gradient of the I/z(t) curves still show 
very good agreement. 

The agreement is excellent over the range of orbits 
tested, which validates that NBODY6df can reproduce the 
expected angular momentum loss at different radii for both 


5.2 Comparison with GADGET 

5.2.1 Cuspy Models 


The rapid relaxation of low-A models means that even if 
the dynamical friction force exerted on the cluster is correct 
at any epoch (as shown in section 5.1.21, different cluster 
realisations will diverge in agreement due to the stochastic- 
ity of the relaxation process. As such, one would ideally like 
to perform NBODY6 runs with higher particle number to 
reduce this effect, but the computational cost is too high at 
the time of writing. Simulation nblk took 7 days to run on 4 
GeForce GTX 780 GPUs and 16 CPU cores. Increasing the 
particle number by a factor of 10 would take over 2 years to 
compute. As an alternative we used the softened tree-code 
GADGET to simulate a larger particle number. We would 
like to stress that in agreement with Kim & Morris (20031, 
we cannot accurately describe the internal dynamics of these 
clusters with GADGET. However tidal stripping still occurs 
in a natural way, and we can compare the bulk properties 
(i.e. mass and position) in order to test the validity of our 
dynamical friction perturbation. 

Softened simulations are expected to exhibit slower dy¬ 
namical friction than collisional simulations. Whilst soften¬ 
ing helps with numerical stability and computational speed. 












































10 J. A. Petts, A. Gualandris, J. I. Read 



Figure 12. Distance of the cluster with respect to the Galac¬ 
tic Centre as a function of time for NBODY6df (blue line) and 
GADGET (green dashed line) simulations. 



Figure 13. Distance of the cluster with respect to the Galactic 
Centre as a function of time for dfalOe (blue line) and gtalOe 
(green dashed line). 


it suppresses the close interactions required for dynamical 
friction. It is true that the GADGET simulations have a 
greater mass resolution for the background than nblk, yet 
the eccentricity of the orbit grows faster due to numerical in¬ 
accuracies in the integration. GADGET’S integrator is only 
accurate to 2nd order, and so accumulates the errors dis¬ 
cussed in Section ISTT] faster than NBODYb’s integrator, 
which is accurate to 4th order. The tree’s force calculation 
is also not as accurate as direct summation, with smooth¬ 
ing effects causing the potential to deviate from spherical 
symmetry. This deviation is evident in Eig. |12| which shows 
the inspiral of cluster models at 10, 5 and 2.5 parsecs in 
NBODYhdf and GADGET. The correspondence is much 
better between simulations dfa2.5 and gta2.5 because the 
number density of the background at 2.5 pc is approximately 
30 times greater than at 10 pc, so the simulation effectively 
has a higher resolution background, and appears more spher¬ 
ical when calculated with the tree. With this taken into ac¬ 
count, the agreement between the NBODYBdf and GAD¬ 
GET is rather good. To see if the discrepancy is due to 
low-Y, we ran gta5 with twice the mass resolution of the 
background. The results were indistinguishable and we can 
be convinced we converged on the solution. 

Fig. shows simulations with the same initial con¬ 
ditions as dfalO and gtalO in Fig. |12[ but with an initial 
eccentricity of 0.3. The same general trend is seen as com¬ 
pared with the circular case, in which GADGET inspirals 
slower, but the first few orbits give very good agreement. 


5.2.2 Shallow Models 


In section we discussed how Chandrasekhar’s dynamical 
friction formula has proven to fail in shallow cusps (i.e. when 
7 approaches 0). The reason for this failure is due to the 
assumption that dynamical friction is a local process, thus 
when p{r) is constant, no special treatment occurs. |Read| 
et al. (20061; Inoue (20091 and Goerdt et al. (20101 show 


that when orbiting within the scale radius of a shallow pro- 
Hle, satellites experience an enhancement of the dynamical 
friction force, followed by an abrupt stalling near the core. 


The enhanced force is thought to occur due to super¬ 
resonance with orbits throughout the core (for the link be¬ 
tween Chandrasekhar’s formula and the resonant nature 
of dynamical friction see also [Tremaine fc Weinberg|p~984| 


Weinberg 19861. This causes the satellite to interact with 


each background star more than once, leading to greater dy¬ 
namical friction as compared with Chandrasekhar’s deriva¬ 
tion, in which the satellite interacts with each background 
star only once. It should again be noted that this is why 
Chandrasekhar’s formula works so well for cuspy profiles, 
where the density is a strong function of R^. Whilst migrat¬ 
ing to the centre of a cuspy profile the satellite constantly 
experiences a whole new set of resonances. These new res¬ 
onances act independently of previous interactions, and as 
such interacting with any star only once is a reasonable ap¬ 
proximation. 

Read et al. ( 2006|) showed that Chandrasekhar’s for¬ 


mula with a fixed log(A) fails at reproducing inspiral in N- 
body models of a harmonic core. However it may be possi¬ 
ble to reproduce these effects with a varying Coulomb loga¬ 
rithm. 

In Fig. we compare the 7 = 0 simulation dfgO.O with 
a self consistent GADGET simulation gtgO.O. We choose to 
only compare GADGET simulations for the shallow models 
because we cannot compare the full inspiral in NBODY6 due 
to the fast relaxation effects at low N, as seen in previous 
sections. Studying the full inspiral for the shallow proHles 
is important as they show interesting deviations from the 
standard Chandrasekhar’s formula. 

Fig. shows that this prescription cannot fully re¬ 
produce the full inspiral of the cluster during the “super- 


Chandraskehar” phase (Read et al. 2006 Goerdt et al. 


20101. The nature of this enhanced force is discussed in sec¬ 


tion Interestingly the stalling radius is very well repro¬ 
duced, thus the physics of the stalling is well captured in 
this model. 

The cluster stalls rapidly due to a combination of two 
correlated effects. Firstly, the low circular velocity in the 
core means that 6min is dominated by GM/v^typ as opposed 
to the cluster size and becomes very large and comparable 























































A semi-analytic dynamical friction model that reproduces core-stalling 11 




Figure 14. Distance of the cluster with respect to the Galactic 
Centre as a function of time for dfgO.O (blue line) and gtgO.O 
(green line). 


to fomax- The physical interpretation in Chandrasekhar’s for¬ 
malism is that any background particles passing by the clus¬ 
ter are deflected by more than 90 degrees, and thus stochas¬ 
tic changes in orbit due to two-body relaxation dominate 
over the many-body dynamical friction effect, suppressing 
any further inspiral. These stochastic changes in orbit are 
evident in the full A'^-body model. The black dashed line in 
Fig. |15| shows the minimum impact parameter fo r the 7 = 0 
case if one instead defines Utyp^ = 2 (j^ -l-u^s as in 


Just et al. 


(20111. It can clearly be seen that the stalling effect is not 
correctly captured if this prescription is used, as the impact 
parameters are not equal until Rg ~ 0. In cuspy profiles 
most of the mass is contained within a small radius, and 
stalling occurs very close to the centre of the system. 

Secondly, as the velocity dispersion doesn’t decrease as 
rapidly as the circular velocity, and is non-zero at J?g = 0 , 
X = vs/V2a shrinks as Rg goes to zero. As such the frac¬ 
tion of stars moving slower than the satellite (those that 
contribute to dynamical friction) decreases rapidly as the 
cluster approaches the core. Fig. shows the fraction of 
stars moving slower than the circular velocity as a function 
of radius for 7 = 0 and 1.5. In the 7 = 0 case very few 
stars move slower than the circular velocity in the core, and 
thus very few stars contribute to dynamical friction. In the 
cuspy case this fraction of stars does not go to zero, and 
the stalling effect is dominated by log(A) approaching zero, 
which occurs at very small radii. If the orbit is initially ec¬ 
centric, faster velocities are reached at pericentre, decreasing 
femin and causing the satellite to stall closer to the centre of 
the host, as can be seen in fig. 3 of Read et al. (20061. 

Fig. shows a similar setup for the 7 = 0.5 profile, 
which shows excellent agreement with the A-body simula¬ 
tion. A cusp of 7 = 0.5 is enough to keep /(u*) from going 
to zero, and the stalling is again dominated by log(A) —> 0. 

It is thus evident that cluster stalling at large radii is 
characteristic of profiles that contain a constant density core, 
such that the cluster stalls when: 


( 10 ) 


Figure 15. Maximum and minimum impact parameters as a 
function of radius for 7 = 0 ,1.5, up to the scale radius. It should 
be noted that the minimum impact parameter is a function of 
the mass of the satellite and background, here the values used 
are the same as throughout this paper. The dashed black line 
shows the minimum impact parameter if the prescription from 
[Just fc Penarrubia] l |2005[ | is used for 7 = 0 . 



Figure 16. Fraction of stars moving slower than the circular 
velocity as a function of radius for 7 = 0,1.5. The distribution 
of stellar velocities is assumed to be Maxwellian with velocity 
dispersion cr, as in equation^ 


and dynamical friction ceases when either Ms ^ Mg{Rg) 


or Vc 


(when 


O.Sa, /(iJ*) ~ 0.03). As both terms 


decrease with decreasing satellite velocity, in profiles with a 
true isothermal core, satellites can stall at Ms/Mg{Rg) ^ 1 
(see section 6 . 11 . 


5.3 N-dependence study with NBODY6df 

We ran a series of simulations to see if the expected N- 
dependence of cluster relaxation, and its effects on inspi¬ 
ral, are well reproduced by NBODYGdf. Simulations dflk- 
df64k have the same initial cluster mass split evenly amongst 


log(A)/(w*) -s- 0, 
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Figure 17. Distance of the cluster with respect to the Galactic 
Centre as a function of time for dfgO.5 (blue line) and gtgO.5 
(green line). 


their cluster particles, and otherwise have the same initial 
conditions. Low-A'^ systems should lose their mass faster 
than an equivalent realistic cluster due to shorter relaxation 
timescales. This behaviour is illustrated in Fig. |18[ which 
shows the bound mass as a function of time for simulations 
with different N. 

With 32k particles the cluster initially has a relaxation 
timescale of ~ 2Myr, and the mass loss is mostly dominated 
by the tides during inspiral (Fig. 18 1 . In simulations df32k 
and df64k, most of the mass is lost when the cluster reaches 
the centre of the potential, where the remaining mass is 
deposited in a disk, and is formally bound to the cusp as 
opposed to itself. 

If the cluster was modelled in a realistic fashion with 
a mean mass of 0.58 Mq (Kroupa 20011 the relaxation 
timescale would be longer than the inspiral time (T^eiax ~ 
9Myr), and the mass loss would be dominated by the shrink¬ 
ing tidal radius. The low-Af models show accelerated mass 
loss due to increased dynamical ejections as expected. 

Fig. |19| shows how this mass loss drastically alters the 
evolution of the orbit. If the cluster loses significant mass, its 
inspiral will stall due to a continually decreasing dynamical 
friction coefficient. 


Figure 18. Mass enclosed in the Roche volume as a function of 
time for simulations dflk-df32k 



Figure 19. Distance of the cluster with respect to the Galactic 
Gentre as a function of time for simulations dflk-df32k. 


6 DISCUSSION 

In this section we consider the inspiral of a point mass ob¬ 
ject, and as such take our dynamical friction formalism and 
implement it as an external force in an independent 2 nd- 
order integrator, which integrates the motion of a point par¬ 
ticle in a Dehnen potential. 


6.1 Understanding stalling in cored profiles 


Read et al. (2006 \ studied the stalling of models that contain 


intrinsic cores using the alpha-beta-gamma prohle (Zhao 


19961: 


p{r) 


Po 

{r/ap{l + (r/a)“j('3-7)/Q ’ 


( 11 ) 


with Po — 9.93 X 10^M 0 kpc“^, a = 0.91 kpc, a = 1.5, 
/3 = 3.0 and 7 = 0.0. The authors could not repro¬ 
duce the stalling effect with a semi-analytical form of 
Chandrasekhar’s formula, as their model used a constant 







































A semi-analytic dynamical friction model that reproduces core-stalling 13 



Figure 20. Semi-analytical integration of 2 X 10® Mq point-mass 
cluster in a uniform background density distribution with param¬ 
eters similar to [Read et al.| l|2006[l . 


Coulomb logarithm. In their model the cluster still would 
have stalled if they had integrated for long enough, due to 
f{vt) —> 0. However this would have happened at much 
smaller radii, as they use a constant Coulomb logarithm. 

If we consider a toy model where we approximate the 
inner region of the core as a constant sharply truncated at 
a, the Jeans equation for the velocity dispersion becomes 
i Binney &: Tremaine|2008 1 : 


vf{r) 



( 12 a) 


where vlr) is the number density and is the potential. We 
can then write: 
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and: 


X = 


\J‘l(o? — r^) 


( 12 g) 


Therefore, Vc < Vr for r < aj^/2, leading to a quickly van¬ 
ishing f{vt) when approaching the core, causing the cluster 
to stall much further out. This coupled with the shrinking 
Coulomb logarithm causes stalling behaviour unique to the 
core model, in which Mg(f?g) Ms. In Fig. |20| we mod¬ 
elled a 2 X 10® Mq point mass cluster in a constant density 


background with po and a the same as in Read et al. 1 20061. 
The cluster stalls at Ms/Mg(i?g) ~ 0.06 in agreement with 
the Wbody simulations of Read et al. (20061. At ~ 200 pc. 


/(u*) ~ 3x 10“®, i.e. only 0.3% of stars move slower than the 
circular velocity, and dynamical friction practically ceases. 
Their cluster stalls a little further in as their model follows 
the distribution of equation [m and thus a is slightly lower 
at ~ 200 pc. 

Our model relies on the key ansatz that dynamical fric¬ 
tion ceases when bmin ^ femax. This can be understood as a 
lack of particles available to be scattered by less than 90 de¬ 
grees by the satellite, but the true structure of such a state 
is not immediately clear. Read et al. (2006 1 argued that the 


system finds itself in a state in which the background parti¬ 
cles orbit in the combined potential of the satellite and the 
harmonic background, in which the time averaged drag force 
on the satellite can be shown to be zero. The authors showed 
that this state gives an excellent match to the phase space 
structure of background stars in the A^-body simulations. 
Thus we can understand this configuration being exactly 
what happens when femin ^ &max, and the system remains 
in a state where no net momentum is transferred from the 
background to the infalling body. 

The ’’Super-Chandrasekhar” inspiral phase prior to the 
rapid stalling is currently not captured by our model, how¬ 
ever it could possibly arise from the contribution to dynam¬ 
ical friction of the stars moving faster than the satellite. |An-| 
tonini & Merritt (2012 1 showed that in shallow cusps around 


black holes, where most stars are moving faster than the cir¬ 
cular velocity, the contribution from the stars moving faster 
than the satellite is significant. This mechanism is likely to 
be similar for the case at hand, however it could well be the 
case that this contribution is not enough to explain the en¬ 
hanced frictional force, and may result from super-resonance 
inside the core | Tremaine fc Weinberg|lSiM Weinberg|1986 


Read et al. 2006 Goerdt et al. 20101. The contribution of 


the faster stars and their effect on the stalling radius is not 
trivial, and shall be left to future work. 


6.2 Comparison with Arca-Sedda &: 




Capuzzo-Dolcetta ^2014 ) 

|Arca-Sedda Capuzzo-Dolcetta| ( |2014[ ) (hereafter AC 14) 
studied dynamical friction in cuspy galaxies and presented 
a new treatment for massive objects near the centre of 
their host systems. The authors derive a semi-analytical for¬ 
mula for the inspiral time of massive point particles orbiting 
Dehnen models, calibrated by A-body models in the GPU- 
parallel direct A-body code HiGPUs ( |Capuzzo-Dolcetta| 
et al.|[20l3 1. In their semi-analytic fitting process, they use 


an exponential interpolation between Chandrasekhar’s for¬ 
mula with a constant 6 max and varying 6 min, and their de¬ 
tailed evaluation of the frictional force near the centre of 
host systems. The authors do not fix bmin and instead let it 
be a fitting parameter, along with Vcr, which they define as 
the critical radius at which they switch to the new regime. 

In our model both femax and femin vary along the orbit as 
a function of the local background and satellite properties. 


Arca-Sedda & Capuzzo-Dolcetta (20141 state that the local 


approximation overestimates the effects of dynamical fric¬ 
tion in the innermost cuspy regions of galaxies, however in 
our approach the maximum impact parameter tends to zero 
at small radii for cuspy distributions, reducing the range at 
which the local approximation acts over. For distributions 
with 7 > 1 this local scale length is smaller than the distance 
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Simulation 

Code 

(pc) 

e 

7 

gtptl.5 

GADGET 

5.0 

0.0 

1.5 

gtptl.Se 

GADGET 

5.0 

0.3 

1.5 

gtptl.O 

GADGET 

5.0 

0.0 

1.0 

gtptl.Oe 

GADGET 

5.0 

0.3 

1.0 

saptl.5 

Semi-Analytic Integrator 

5.0 

0.0 

1.5 

saptl.5e 

Semi-Analytic Integrator 

5.0 

0.3 

1.5 

saptl.O 

Semi-Analytic Integrator 

5.0 

0.0 

1.0 

saptl.Oe 

Semi-Analytic Integrator 

5.0 

0.3 

1.0 


Table 2. Simulations in which the cluster is modelled by a point 
mass particle. Simulations are performed in a self consistent way 
in GADGET, the mass resolution of the GADGET simulations 
is 30 M 0 . The properties of the cluster and background are the 
same as in section 



Figure 21. Distance of the cluster with respect to the Galactic 
Centre as a function of time for simulations gtptl.5, gtptl.Se, 
saptl.5 and saptl.5e. 


to the centre of the background. By using 6 max = Rg/'y we 
ensure that the local approximation is valid, but probably 
slightly underestimate the frictional force at the very centre 


of the systems, as AC14 suggest (see also Just & Pefiarrubia 


2005). 


|Arca-Sedda fc Capuzzo-Dolcetta] ( |2014| model the 
satellites as Plummer-softened point particles. Simulations 
with a satellite consisting of a cluster of particles take longer 
to reach the centre of their host, due to mass loss and the 
larger bmin of extended objects. For this reason we cannot 
directly test NBODYGdf against their timescale formula. In¬ 
stead we use our semi-analytical integrator. We also perform 
GADGET simulations with the same initial conditions to 
compare our results. The list of simulations is presented in 


Tabled 


We compared the inspiral time of these simulations 
with results from AC 14 and found significant discrepancy 
with their dynamical friction timescale, which was calibrated 
mostly on radial orbits (equation 20 in AG 14). However a 
good agreement is found with an improved formula cali¬ 
brated on a wider range of models, given in|Arca-Sedda et al.| 
(arXiv:1501.04567). 

Both our semi-analytic approach and GADGET sim¬ 



Figure 22. Distance of the cluster with respect to the Galactic 
Gentre as a function of time for simulations gtptl.O, gtptl.Oe, 
saptl.O and saptl.Oe. 


ulations show good agreement with the revised timescale 
formula, (see Fig. 22 and 21 for the 7 = 1.0 and 7 = 1.5 


cases, respectively). The radial trajectory of the inspiral in 
GADGET is very well reproduced by our semi-analytic for¬ 
mula and we can validate our approach for the inner cuspy 
regions. In the 7 = 1.0 the semi-analytic approach diverges 
slightly from the GADGET simulation, as the live back¬ 
ground distribution is slightly shallowed by the inspiraling 
body, however the match is still reasonably good, with the 
inspiral time being well captured. A mechanism for feeding 
the energy lost by the satellite back into the analytic back¬ 
ground distribution would be able to correct for this effect 
for massive satellites in shallow background distributions, 
however this is beyond the scope of this work. 

It should be noted that gtptl.5e agrees much bet¬ 
ter with our semi-analytic model than gtalOe does with 
NBODYGdf. This is because the effect of close encounters 
(i.e. b ~ bmin) is completely captured by the high mass ratio 
in the point mass case. Within a particular cluster, the effect 
of close encounters with background stars at the edge of the 
cluster are underestimated. The NBODY 6 simulations treat 
these encounters properly, and thus excellent agreement is 
found. 

Chandrasekhar’s local approximation is inaccurate near 
the centre of cuspy host systems, and AC14’s approach is 
more representative of the true force in the very central re¬ 
gion of the background distribution. We recommend the use 
of AC14’s numerically calculated Coulomb logarithm when 
a very accurate representation of inspiral is required in the 
very inner region of cuspy profiles. However, our varying 
Coulomb logarithm can reasonably approximate the force 
experienced throughout the cuspy region, only slightly un¬ 
derestimating the inspiral time, without the need for two 
free parameters. 


7 CONCLUSIONS 

We present a modification to the GPU-enabled direct sum¬ 
mation code NBODY 6 , which we name NBODYGdf, to in- 
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elude the effects of dynamical friction on the inspiral of a 
star cluster in a smooth background particle distribution. In 
this approach, the dynamical friction force on each cluster 
particle is computed analytically and added to the A'^-body 
forces exerted by the other cluster particles. In this way, 
only the cluster needs to be modelled in a A^-body fash¬ 
ion, while the effect of the background stars is included in 
an approximated but reliable way. This significantly reduces 
computational time with respect to a full A'^-body modelling 
of the cluster and the background system. 

It should be emphasised that the dynamical friction 
treatment implemented in NBODY6df is physically moti¬ 
vated rather than calibrated on A^-body simulations, and 
thus has predictive power, owing to the physically motivated 
maximum and minimum impact parameters. The predictive 
power of NBODYGdf allows for quick modelling of a large pa¬ 
rameter space of initial conditions without prior calibration. 
The mass term in Chandrasekhar’s formula for extended ob¬ 
jects is found to be well represented by the mass enclosed 
in the Roche volume, as opposed to just the formally bound 
stars. This is due to the presence of potential escapers en¬ 
hancing the gravitational wake whilst they are still close 
to the cluster. NBODYGdf can be used to simulate young 
cluster inspiral in the Galactic Centre, or the inspiral of 
globular clusters or dwarf galaxies in the halo of a larger 
host. It should be noted that dynamical friction in a disk 
or other highly non-spherical systems cannot yet be reliably 
modelled with NBODY6df. This is due to the maximum im¬ 
pact parameter being smaller perpendicular to the disk than 
parallel to it. Accurately modelling inspiral in a disk would 
require an angular dependence in the summation of possible 
impact parameters, and is beyond the scope of this work. A 
comprehensive study of young dense clusters formed close to 
the Galactic Centre will be presented in an upcoming paper. 

NBODY6df can also accurately model the inspiral of 
satellites in shallow profiles, due to a new approach in which 
the minimum parameter is defined to be inversely propor¬ 
tional to the square of the satellite’s velocity, and indepen¬ 
dent of the velocity dispersion. 

For a direct summation code, computational time scales 
with N'^ for an integration of one A^-body time unit. A full 
A-body simulation of a 10® Mq cluster with mean mass of 
0.58 Mq and a Kroupa mass function would require ~ 9 x 
10^ background particles for a 10:1 ratio of Mbg : Me. On 
the other hand, a simulation with NBODYOdf would only 
require the 1.73 x 10® cluster particles to be modelled as A- 
body particles, reducing the computational time by several 
orders of magnitude. 

While the current implementation adopts a Dehnen 
model for the background system, any static model can be 
implemented in order to follow the evolution of star clusters 
in which dynamical friction of the orbit is important. The 
code will be released publicly on githulQ 
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APPENDIX A: VELOCITY DISPERSION 


Here we derive the velocity dispersion as a function of radius 
for Dehnen models. We also include the optional potential 
of a central black whole, which although not used in this 
paper, is available in the initial release of NBODY 6 df. 

The velocity dispersion for an isotropic spherical system 
is ( Binney fc Tremaine||200^ |: 


Where v{r) is the number density and 4? is the potential. 
The number density of a Dehnen model is given by: 


^ = (3-7) a 

Mg 47r rT'(r + a)'^— 


And the potential by | Dehnen|1993 1: 


<I>(r) = — 


GMg 


a 2 — 7 


1 - 


r 2-7 


r + a 


GfiMg 


(A2) 


(A3) 


Where the second term is the additional potential due to an 
optional central black hole, and is the ratio of the black 
hole and Dehnen model masses. It follows that: 


dr 


GM, 


(*) 


1-7 


(o + r )2 
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G/rMg 


(A4) 


By putting equations |A2| and |A4| into equation | A1 1 and mak¬ 
ing the substitution x = r/a, equation |A1| becomes: 

vA(x) = x^(x + [f{x) + fJ,h{x)] (A5) 


Where: 




x'^{x + 1)®-T' 


dx 


(A 6 ) 


POO 

h{x) = / x~"'~^{x dx (A7) 

J X 

Which must be evaluated for the desired value of 7 . f{x) 
and h{x) are analytic for integer values of 47 . If fj, is zero 
the h{x) term is skipped in NBODY 6 df. 














